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Electronic structure calculations for compressed molecular 
hydrogen are performed to provide more insight into the di- 
versity of phenomena recently observed experimentally. We 
perform full-potential LAPW calculations and analyze them 
in terms of a molecular tight-binding model. We show that 
a g and a u bands overlap occurs at rather low pressure, while 
an insulating state persists at specific orientations up to 200 
GPa, in accord with previous work, and is due to opening of 
hybridization gaps at the Fermi level. We also present cal- 
culations of electronic properties such as plasma frequencies 
and electron-vibron coupling. 
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Experimental data on the behavior of hydrogen at 
megabar pressures is rapidly emerging as a result of ad- 
vances in in diamond-anvil cell tcchique. Certain results 
are clearly established and understood, while a number 
of questions remain, most of which concern the implica- 
tion of spectroscopic results for the electronic structure. 
One issue, in particular, is band overlap metallization. 
It is well established that at room temperature and be- 
low hydrogen under high pressure stays insulating up to 
at least 150 GPaEI (which is ahout. tenfold compression, 
r s sw 1.4). LDA calculationsBBu'El give a metallization 
pressure that strongly depends on the assumed orienta- 
tion of hydrogen molecules (reported numbers are from 
40 to 150 GPa), the strongest metallization tendency be- 
ing in the structure where all hydrogen molecules are 
parallel. The reason for such a strong orientational de- 
pendence is not completely understood. 

A very basic argument for band-overlap metallization 
of hydrogen is that the insulator-metal transition occurs 
when the bandwidths of the lo~ g and la u bands originat- 
ing from the corresponding molecular orbitals become 
larger than the splitting of the two molecular levels. One 
can easily estimate when this should happen using, for 
instance, Slater's variational LCAOEl, and assuming the 
perpendicular orientaticua (which minimizes the overlap) 
of neighboring moleculeaJ. The answer is r s sw 1.45. More 
accurate numerical calculations, starting from Friedli and 
Ashcroft,H, have yielded similar numbers. However, it 
turns out that even LDA calculations (see, e.g., Refs. 
§i (§ |5|) render an insulating ground state in the ori- 
entation which minimizes the total energy. In the local 
density approximation (LDA), where molecular orbitals 
are too diffuse, this sort of metallization occurs at even 
larger r s . 



A detailed microscopic analysis of compressed hy- 
drogen, even in LDA and for clamped nuclei, has not 
been yet presented. To provide further insight into the 
available experimental data, we performed full-potential 
LAPW calculations for hydrogen for different orienta- 
tions. By mapping the calculated band structure onto 
the two molecular orbitals nearest-neighbor tight bind- 
ing model we show how this band structure is formed 
and why the minimum total energy orientation appears 
insulating. n |-||-| 

Most previous calculationsatm used a plane wave ex- 
pansion for the wave functions. This is a good basis set 
for the interstitial region; however, close to nuclei it con- 
verges slowly. For a standard cut-off of 50 Ry the Fourier 
expansion of hydrogen Is wave function converges at a 
nucleus only within 14%. We used the Linearized Aug- 
mented Plane Wave method with an MT radius of 0.69 
bohr; for r > 0.69 the error in the aforementioned Is 
Fourier expansion is about 50 times smaller than for 
r = 0. We found the total energy to be convergent better 
than 1 mRy, with the cut-off of 44 Ry. We find that the 
calculations converge relatively slowly with the number 
of k-points in the Brillouin zone; we had to use a special 
k-point mesh 12x8 x 8 corresponding to 96 points in the 
irreducible edge for the 4-molecule orthorombic structure 
described below. Most calculations reported here were 
performed for r s = 1.435 (lattice parameter 3.3 bohr, 
molecular volume 2.209 cm /mol, p/po — 10.5). The cal- 
culated pressure was about 150 GPa at this compression. 
We also performed some calculations at higher compres- 
sions in order to detect the LDA metallization. 

The crystal structure of high-pressure phases II and 
III has not been determined^. It is generally believed 
that centers of molecules form an hep lattice, as they are 
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known to do to at least 50GPa. At least four distinct 
optical (Raman and infrared) vibron mades are observed 
in phase II in hydrogen and deuteriumElfl, indicating an 
orientational ordering of at least four molecules per cell. 
Phase III has an orientational order different from that 
of phase IIO, with onepStrong Raman and one extremely 
strong infrared vibronlH. Both vibrons show substantial 
softening, typically by 200 cm -1 between 200 and 100 K, 
when passing into phase III. The structure considered the 
most relevant for the high-pressure phases of hydrogen is 
the Pca2i structure which is derived from the hcp-c (in 
which all molecules are aligned parallel to c, space group 
P63 / mine) by doubling the unit cell and making the four 
molecules in the unit cell to be tilted from c by an angle 
9 and rotated by ±ip, it ± (p. This structure is known to 
minimize the electrostatic energy of classical quadrupoles 
on an hep lattice, at 9 ss 55° and ip w 43°. Notably, even 
at p/po ~ 13 (r s sa 1.35) the equilibrium orientation, 
calculated from total energy methods, is similar to thisu, 
although quadrupolar interaction is is not the dominant 
anisotropic interaction at this compression. Our compu- 
tations, described below, give 9 ~ 62° and tp w 50° at 
r s = 1.435. 

The assumed structure has five degrees of freedom: 
c/a and b/a, two orientational angles, and the H2 bond 
length. We have fixed b/a = as in hep, and optimized 
the structure with respect to four remaining parameters. 
We find that c/a is not sensitive to orientation: for the 
high-symmetry PGa/mmc structure, i.e. for 9 = 0, we 
find c/a = 1.59, while for (9,p) = (62°, 50°) orientation 
it is 1.62. The optimal bond length changes more: for the 
same c/a (1.59) the bond length d in P63/mmc structure 
is 1.42 bohr, while for (62°, 50°) it is 1.48 bohr. Most of 
our calculations below use c/a — 1.6 and d = 1.5 bohr. 

Fig.|l| shows the orientational dependence of the to- 
tal energy. Note that the potential barrier for librations 
around (62°, 90°) is rather low, about 280 K. Of course, 
the calculated energy profile corresponds to coherent ro- 
tations. Actual behavior of the system depends on how 
the energy profile differs when a single molecule is rotat- 
ing in the field of its neighbors. For a single molecule 
librating in the potential shown in Fig.[l], we can solve 
the Schrodinger equation for the anharmonic oscillator 
numerically and find that even at zero temperature the 
molecule would be freely librating around the average 
orientation (62°, 90°) with the amplitude of «40° (both 
for hydrogen and deuterium), thus making a less orien- 
tationally ordered phase. If we want to characterize the 
degree of orientational order, we should introduce an ori- 
entational order parameter, which is usually written as 
a zero-trace tensor Q a p — (3(n a np) — 5 a p)/2, where n 
is the unit vector of the molecule's direction. Principal 
values of this tensor for axially symmetric librations are 
{Q,-Q/2,-Q/2}, and Q = 1 for full ordering. Libra- 
tions described above correspond to Q aa « {0.5, 0, —0.5}, 
while free rotation around c-axis with 9 = 62° gives 
Qaa = {-0.17,0.085,0.085}. 



The first two azimutal librational eigenstates are about 
30 cm -1 (7.5 cm -1 ) and 170 cm _1 (100 cm -1 ) above the 
saddle point at (62°, 90°) for hydrogen (deuterium). 
Since (62°, 90°) is a saddle point, azimutal and po- 
lar librations are quite different: Azimutal librations, as 
described above, are strongly anharmonic, while polar 
vibrations, to the contrary, can be treated in the har- 
monic approximation. Their amplitude reaches ~ 40° at 
T ~ 400K. The calculated frequency of the polar libron 
is about 350 cm -1 . 

One of the most intriguing issue is whether or not the 
equilibrium geometry is governed by electrostatic forces. 
The calculated ground state (62°, 50°) is very close to the 
equilibrium orientation (55°, 43.5°) for the quadrupolar 
hep lattice (see, also, Ref. [|). This orientation (62°, 50°) 
maximizes the insulating gap (actually even small de- 
viations from this orientation render a metallic state). 
Thus, the band structure energy also favors this orienta- 
tion. Besides, at the tenfold compression the actual elec- 
trostatic interaction differs substantially from the ideal 
quadrupolar one. These arguments cast doubt about 
quadrupolar interaction being the main ordering factor. 
More insight can be gained from comparing calculations 
with different number of k-points. We have performed 
calculations with a sparser mesh (12 points) and found 
that the total energy is minimal at (20°, 90°). 12 k-points 
for an 8-atoms unit cell is a sufficiently large number to 
have the charge density converged. Indeed, we found no 
substantial difference between the orientational depen- 
dence of the electrostatic energy in the both sets of calcu- 
lations. To the contrary, the one-electron energy orienta- 
tional dependence is quite different in the two cases (with 
12 points, it is not converged yet). Thus, we conclude 
that the main factor controlling ithe orientational order- 
ing at about tenfold compression is one-electron (band) 
energies. 

Since the driving force for orientational ordering is 
the band energy, it is useful to compare the electronic 
bands for several different orientation. Fig.||a shows the 
bands (full lines) for the two-molecules-per-cell P63/mmc 
hep structure. To understand the nature of the rele- 
vant bands, we compare the self-consistent band struc- 
ture with that from two-molecular-orbital tight-binding 
(TB) model (dashed lines). The TB parameters were 
t = 3 eV, E au - E ag = 16.8 eV at r s = 1.435. Only 12 
nearest neighbors were taken into account. The lowest 
states at T and X are pure a g [actually, cr g (l) ± a g (2)]. 
However, the next level at T is already the bottom of the 
a u band. Thus the overlap of the a g and a u bands is very 
large, of the order of one Rydberg. Note also the third 
band at T, which is the bottom of the 2a g band. 

Since the a g — a u band overlap occurs at much lower 
compression (p/po ~ 7), H2 would be metallic if not for 
the gaps that open at the Fermi level due to the orien- 
tational tilting. Fig.^D shows bands for a less symmet- 
ric (62°, 90°) structure (dashed lines), compared to the 
folded down hep bands (dotted lines). At several points 
where two downfolded bands cross in the hep structure, 



2 



the degeneracy is lifted now because of lower symmetry. 
This band structure is still metallic; the crossings occur 
at different energies. Nevertheless, the density of states 
at the Fermi level and correspondingly the one-electron 
energy is decreased. 

If we now move from the (62°, 90°) structure further 
to the equilibrium orientation (62°, 50°), (full lines in 
Fig. ||b) we observe that local minima in the conduction 
band, and maxima of the valence band (which emerged 
from the band crossings in the P63/mmc structure) are 
now aligned at, respectively, 18 eV and 17 eV, thus form- 
ing a considerable gap. Rotating the H2 molecule influ- 
ences the positions of these maxima and minima in a 
complicated way; alignment of the direct gaps such that 
an indirect gap appears is characteristic for the equilib- 
rium orientation. There is an analogy with the Peierls 
transition in metals with Fermi surface nesting: if a su- 
perstructure exists in which the nesting bands become 
degenerate, lowering of symmetry removes the degen- 
eracy, decreases the band energy and induces a struc- 
tural phase transition. An important difference is that 
there is no nesting, and therefore the regular suscepti- 
bility would diverge in the high-symmetry phase, as in a 
Peierls transition. Instead, there is an instability with re- 
spect to final rotation of the molecules. Correspondingly, 
the known problem of the one-dimensional Peierls tran- 
sition, namely that the divergence in the non-interacting 
susceptibility may cancel out in the full susceptibilitytil, 
does not apply. 

Another way to rationalize this result is to use Ander- 
sen's force theoremO, which states that the first-order 
total energy differences in LDA can be calculated as a 
sum of the electrostatic energy differences and band en- 
ergy differences, provided that charge density is not re- 
calculated to self consistency in the perturbed configura- 
tion (for instance, if one would ascribe a sphere around 
each molecule and rotate it rigidly without changing the 
charge distribution inside the sphere). While it was not 
done in our calculations (and is technically not possible 
in LAPW method), the changes of molecular charge den- 
sity upon rotation are small. 

We have also performed calculations for higher com- 
pression, but without orientational optimization. The 
indirect gap closes after further .compression by 6.5% 
(P w 185 GPa). It is well knowntj that the LDA con- 
siderably underestimates semiconductor gaps. Therefore, 
not only should hydrogen stay insulating in the pressure 
range of 150 GPa, but the energy of the gap is likely to 
be rather large. Of course, dynamical disorder smears 
the gap, thus one can expect gradual increase of the low- 
energy absorption with temperature and pressure in a 
non-Drude manner, with eventual full metallization at 
pressures larger than at least 185 GPa. Optical exper- 
iments do not r . ev .eal metallic, Drude-like absorption up 
to -250 GPaMlla; In the phase III which exists above 
150 GPa, strong activity of the infrared vibron with fre- 
quency i/«4x 10 3 |Sm _1 (in deuterium « 3 x 10 3 cm -1 ) 
is well documented^. This puts an upper limit to the 



plasma frequency: if hydrogen is a semi-metal at this 
pressure, then lo v \ < 0.4 eV. 

The main vibron frequency softens discontinuously 
at the transition to the phase III, which in some pa- 
pers was attributed to the onset of metallic screening I , 
and keeps softening with cooling (up to ~ 200 cm _1 )E. 
Current calculations provide several strong arguments 
against this pointjaLjview. First, as also noticed in pre- 
vious calculationsotffl, an arbitrary orientation results in 
a smaller gap or in the closure of the gap, compared to 
the optimal orientation. Thus one can expect that at ele- 
vated temperature when the molecules do not assume any 
specific orientation, metallization occurs at lower pres- 
sure, and upon heating, not cooling. This would result 
in hardening of the vibrons with cooling, and not soft- 
ening, and would also demand that phase boundary line 
would shift to lower temperature with pressure, opposite 
to the observed trendem. 

It is also instructive to estimate the possible ef- 
fect of metallic screening on vibron frequencies. To 
do that, we can calculate coupling of the main vi- 
bron with conducting electrons, A = —Alu/uj = 
J2nk( de nk/dR) 2 /2Muj 2 S(e nk - E F ), where Auj is the 
change of the vibron frequency u> due to metallic screen- 
ing, R is the H-H bond length, M is the proton 
mass and the deformation potentials (de n k/ dR) are 
the derivatives o£-pne-electron energies with respect to 
the bond lengthllj. We have done such calculations 
for the (<p,0°) and for (65°, 90°) orientations, using 96 
nonequivalent first-principle k-points interpolated onto 
the denser mesh of A90 nonequivalent points, using a 
Fourier interpolation^]. By shifting Ep in the valence 
band we imitate the effect of different degrees of metal- 
lization without doing calculations at the different pres- 
sures. Simultaneously we can control the degree of met- 
allization by calculating the plasma frequency uj^ = 

(8ire 2 /h 2 m e V ce u) X^„k( 5e «k/c'k) 2 (5(e nk - E F ) and com- 
paring it with the above-mentioned upper bound. The 
results are shown on Fig.| One can see that the effect of 
possible metallization on the phonon frequencies is much 
too small to explain the softening. 

To summarize, we presented detailed full-potential 
LAPW calculations for clamped-nuclei molecular hydro- 
gen in the Pca2i crystal structure for pressures of 150 
GPa and above. The structure was optimized with re- 
spect to c/a, H-H bond length, and molecular orienta- 
tions. We found that orientation that minimizes total 
energy is close to the equilibrium orientation of classical 
quadrupoles, although the main driving force for orienta- 
tional ordering is one-electron (band) energy. Molecular 
hydrogen remains insulating well above the compression 
at which a g — a u bands start to overlap, because at the 
equilibrium orientation a g — a u band crossings occur close 
to Fermi energy, so that the hybridizarion open a dielec- 
tric gap. On the other hand, dynamical disorder due to 
zero-point motion will transform the dielectric gap into a 
pseudogap. We also have calculated electron-vibron cou- 
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pling constant and found that in cannot account for the 
5% vibron softeneing in Phase III, which must be related 
primarily with new orienational ordering in this phase. 

We thank A. Goncharov, M. Li, and H.-K. Mao for dis- 
cussions of the high-pressure experiments on hydrogen, 
and particularly R. Hcmlcy for numerous useful com- 
ments. Computations were performed on the Cray Y-MP 
at the NCSA and on the Cray-90 at Pittsburg Supercom- 
puter Center. 



FIG. 2. (a) LAPW energy bands for the hep structure, 
P63/mmc, compared with a (dash lines) two-molecular- or- 
bitals tight-binding model, (b) LAPW energy bands for Pca2i 
structure in (Q°,ip) orientation (downfolded P63/mmc struc- 
ture, dotted lines), (62°, 90°) (dash lines) and (62°, 50°) ori- 
entations. 

FIG. 3. Plasma frequencies and electron-vibron coupling 
constant for (62°, 90°) orientation. Experimental upper limit 
for huip, ?s 0.4 eV, is also shown. 
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FIG. 1. Orientational dependence of the total energy of 
Pca2i hydrogen. Radial coordinate, 6, is the polar angle and 
angular coordinate, ip, is the azimutal angle. Lines are guides 
for eye. 
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